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ABSTRACT 



Four time differencing schemes were tested using a 
barotropic primitive equation model on a spherical 
staggered grid with an analytic input in order to com- 
pare amplitudes, phase speeds, and computation time for 
each. The methods tested were the leapfrog, Euler- 
backward, leapfrog-trapezoidal, and Adams-Bashford. One 
set of experiments was performed using an averaging 
technique to reduce the effects of gravity waves in the 
higher latitudes. Another set was performed without the 
averaging in order to determine the effects of this 
technique on the solutions. 
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- Arbitrary constant in the stream function 

- Arbitrary constants for fourier series cosine terms 
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- Earth's radius 
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- Arbitrary constants for fourier series sine terms 

- Arbitrary constants for fourier series combined terms 

- Wave speed 

- Grid increment 

- Constant used in the Arakawa averaging method 

- Greatest integer value of Dj 

- Arbitrary function 

- Coriolis parameter 

- Acceleration of gravity 

- Height of free surface 

- Height of free surface at wind points 

- Height computed by equation (3) in "restorative- 
iterative" method 

- Height derived from the linear balance equation and 
used as an "observed" value in the "restorative- 
iterative" method 

- Grid index in the x-direction (east-west) 

- Grid index in the y-direction (north-south) 

- Restoration coefficient for zonal wind 

- Restoration coefficient for meridional wind 

- Restoration coefficient for height 

- Variable grid index in the x-direction which is 
dependent on j 
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X - Longitude 

v - Angular wave speed 

ip - Stream function 

V - Del operator (horizontal) 

- Laplacian operator (horizontal) 
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I . INTRODUCTION 



In the field of operational numerical weather prediction, 
the trend, in recent years, has been toward the development 
of sophisticated global prediction models. This has been 
made possible by the rapid expansion of computing capacity 
and developments related to general circulation research. 

The purpose of this study was to examine various time 
differencing methods, using a barotropic primitive equations 
model on a global staggered grid. A spherical harmonic 
analytic stream function was used for the initial conditions. 
By using an analytic initial condition, errors in real data 
observations and analysis, which are unavoidable in practical 
dynamical prediction, were eliminated. 

The objective was to compare the time required for com- 
putation, amplitudes, and phase speeds for each of the time 
differencing schemes. 
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II. BAROTROPIC PRIMITIVE EQUATIONS MODEL 



Time differencing experiments were performed using the 
free surface barotropic primitive equations. The integra- 
tions were carried out on the sphere using the difference 
method of the Arakawa type which was developed by Winninghoff 
(1971) . 

A. PRIMITIVE EQUATIONS 

The primitive equations , in spherical coordinates and in 
flux form, for this model are: 



3 (uh) 



3t 



a cos 0 



3(uuh) 3(uvh cos 0) 

3 A 30 



uvh tan 0 

+ + fvh - 



a (gh) 



a cos 0 3 A 



(1) 



3 ( vh) 



3t 



a cos 0 



3(vuh) 3 (vvh cos 0) 

3 A 30 



uuh tan 0 h 3(gh) 

+ f uh - _ 



a 3 A 



( 2 ) 



3(h) 



3t 



a cos 0 



3(uh) 3(vh cos 0) 

3 A 30 



(3) 



Equations (1) and (2) are, respectively, the zonal and merid- 
ional momentum equations, and equation (3) is the continuity 
equation . 
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B. GRID 



The spatial finite differencing was performed on a 
staggered, spherical grid. The wind and height variables 
were carried at alternate points (see Fig. 1) with height 
only at the poles. The latitudinal and longitudinal grid 
increments were five degrees. This gives 2560 points 
(72 x 35) over the globe, with wind and height carried at 
1260 points each. 

h (North Pole ) 



u,v h u,v h h 

h u,v h u,v u,v 



u,v h u,v h h 

| h u,v h u,v . . • u,v 

j u,v h u,v h h 



i — h (South Pole) 

FIGURE 1. LOCATION OF VARIABLES 

C. SPATIAL FINITE DIFFERENCING 

Spatial differencing of the Arakawa (1966) type was 
used which eliminates the spurious energy growth which can 
occur with standard finite difference approximations to 
the nonlinear advection terms. The difference equations 
used here for this model are: 
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A (uh) ij 



At 



a cos 0j 



(^i j^i+l j ) ( u ^) i+Jjj ( u ij +u i— lj) (uh) i-^j 



2AX 



+ . ^ij^ij+ 2 1 (vh) ij+i cos 9 jti- (u ij tu ij-2) (vh) ij-i c ° s 9 j-i 



2A0 



U^-s Cvh)^ tan 0- 

+ -3 3 1 + fjivh)^ 



gh- 


h i j " h i-lj 


a cos 0j 


AX 



(4) 



A (vh) 



ID 



At 



a cos 0- 



( v ij+ v i+ lj) (uh) i+i-j- ( v ij +v i-lj ) (uh) 



2AX 



<rij +v i j+ 2> < vh >ij+l cos e j + l- (v ij +v ij- 2 ' <vh »ij-lC°s Sj.i 



2A0 

u ii ( uh ) ii tan 0 i 9 h i-i 

— d td L - f . (uh) . . - — il 

D ID 



h l-l j-l _h £-l i-1 



A0 



(5) 



Ah; 



At 



a cos 0- 



(uh) i+ ij - (uhJi.ij 



AX 



(vh) £ j + i cos ©j + l “ ( v h ) £ j — 1 cos ®j-l 



A 0 



( 6 ) 



where 



<uh) ij E u ij h ij 



and 



( vh ) i j = v^h^ 



where 



^^ij + ^i-lj) “ -j-T^l+lj j ** *H-lj + ^i-2j^ 



1 D 



lj 



which is a second order, one-dimensional Bessel's interpola- 
tion scheme with p = h. 



13 



Similarly 



(uh) i + ^ j - h. [ (uh) j+ (uh) j ] - __t (uh) ^ + 2 j- (uh) j 




and 





Originally, all the interpolated, starred quantities were 
derived using a two-dimensional linear interpolation. This 
method was found to introduce a 2d-wave which in certain 
areas was an order of magnitude greater in amplitude than 
the zonal wave used for the experiments. 

Figure 2 shows the indexing convention used for the 
mass and wind variables. The index 1 in equations (5) and 
(6) was i if j was odd and i + 1 if j was even. 



j odd 




j even 





j odd 





FIGURE 2. GRID INDEXING 
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III. TIME DIFFERENCING METHODS 



Four time differencing methods were used to evaluate the 

phase angle and amplitude errors of each method. The errors 

/ 

were evaluated by comparison to an analytic solution to the 
non-divergent barotropic vorticity equation. 

The four methods tested were the leapfrog, Euler-backward, 
leapfrog- trapezoidal , and Adams-Bashford schemes. Two tests 
were made with each method. 

The first set of tests used time increments of fifteen 
minutes for the first three methods and ten minutes for the 
Adams-Bashford scheme. 

These time steps were possible, even though the i grid 
distance at 85° N and S is only 40 km, because of a pro- 
cedure, used by Arakawa (Langlois and Kwok, 1969) to average 
the effects of high frequency inertial gravity waves in the 
zonal direction. 

The averaging technique involved a coefficient 
Aj = . 125 (Dj - 1)/Dj (7) 

where 




and Dj was the greatest integer value of Dj . 

The uh terms in equation (3) and in the advective terms 
of equation (1) and (2) were replaced by 

(uh) 1 = (uh) +A [ (uh) +(uh) . -2(uh) .] 

iD ij j i+lj i-lj ID 

For 1< Dj< 2 (8) 
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and 

(uh) N = (uh)N-l+A [ (uh)N-l.+(uh)N-l .- 2 (uh)N 7 l] 
ij 13 3 1+I3 1-I3 13 

for N<Dj<N+l ( 9 ) 

Similarly h in the pressure gradient term in equation ( 1 ) 
was replaced by 

h 1 = h +A [h +h - 2 h ] for l<D -<2 ( 10 ) 

ij ij j i+lj i-lj ij J 

and 

h 1 ? . = h^-l+A . [h^-l +hN“l -2h^Tl] for N£D-;<N +1 ( 11 ) 

13 13 3 1+I3 1-I3 13 3 

No averaging was done if Dj<l. 

The second set of tests were run without using the 

averaging technique. Thus to remain within the von Neumann 

linear computational stability criterion (Haltiner, 1971 ) , 

a 2 . 5 -minute time step was used. The two sets of tests 

were run in order to determine the effects of the averaging 

technique on the solutions. 

A. LEAPFROG 

The leapfrog method is a centered time differencing 

c A t 

scheme which is conditionally stable for < 1 . The finite 

Ax 

difference equation is: 

F t+1 _ F t-1 + 2 At 9Ft (12) 

9 t 

Since this method has three time levels, it has both a 

# 

physical and a computational mode. 
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B . EULE R- BACKWARD 

The Euler-backward method is a two-step iterative 



scheme which is conditionally stable for <1. The 

Ax 

difference equations are: 

3F fc 

F* = + At 

at 

aF* 

F t+ 1 = Ft + At (13) 

at 



Since the Euler-backward method has just two time levels 
it has only a physical mode. 



C. LEAPFROG-TRAPEZOIDAL 



The leapfrog-trapezoidal is another two-step iterative 

scheme which is conditionally stable for <yJT. 

Ax 



The difference equations are: 



aF fc 

F* = F fci l + 2A t 

at 

At / 3F* 

Ft+1 = F t + — [ + 



2 



at 




(14) 



Since this method, like the leapfrog, has three time levels, 
it also has both a physical and a computational mode. 



D. ADAMS-BASHFORD 

The Adams-Bashford method used was the one examined by 
Lilly (1965) . 

The difference equation is: 

/ 3 3F fc 1 3F t_ l\ 

F fc +1 = F fc + At [ I (15) 

V2 at 2 at J 
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This method has three time levels, thus it has both a com- 
putational and a physical mode. The method is unstable but 
has some desirable features. The computational mode tends 
to damp and the rate of erroneous amplification of the 
physical mode is small if At is small. 
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IV. INITIAL CONDITIONS 



The initial velocity and height fields for these ex- 
periments were derived from a stream function which is a 
solution to the non-divergent barotropic vorticity equa- 
tion. The stream function used was examined by Gates 
(1962) and Neamtan (1946) , which is: 

rp = A sin(mA-vt) P m (sinO) - B a? sin0 + C (sin 0) 

n (16) 

A reasonable meteorological pattern was obtained from equa- 
tion (16) by selecting 
C = 0 

A = 1000 m^ sec - l 

The constant B was related to the angular wave speed by 
V n(n+l)-2 29 

— = B - (17) 

m n(n+l) n(n+l) 

For wave number 6 and, with n = 7 for convenience, — « 20° 

m 

long, per day 

B = 6.8905 x 10“6 sec -l 
The stream function then became 

% 

ip = -279.68 x 10 - 6sinG + 136.65 x 10“^ sin ( 6 A-vt ) sin0 

cos ^0 m^sec"! (18) 

Since these experiments were performed using a free sur- 
face barotropic primitive equations model which allows di- 
vergence, equation (17) was satisfied only approximately. 
Rossby (1939) has shown that the presence of divergence in 
a barotropic atmosphere will slow up the rate of wave propa- 
gation, especially for small values of wave number m. 
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The initial wind field was a non-divergent wind given 



by 



u = 



3i£i 1 9t|> 

9g a 90 



(19) 



9ip 1 9i|> 



v = 



( 20 ) 



9x a cos 0 9A 



The initial height field was derived by solving the 
linear balance equation 



V 2 h = _ [fV 2 i|> + V^*Vf] 

g 



( 21 ) 



where 



■i[ 



1 9 2 1 9 9 

+ — (cos 0 ) 

cos 2 0 9A 2 cos 0 90 90 



] 



and 



v -(■ — 

\a cos 0 9A a 90 / 



Equation (21) was solved by the following relaxation 
scheme : 



h - 1 = 4 . ♦ Rij 



( 22 ) 



where 



.9 



R ii ' 



(COS0-+ ) 

J COS0j 



[ (COS0 . - 
3 



.-d sin0-; ) h . . , _+ (cos0 -+d sin0^)h. . 9 

] J 1J + Z j J -LJ ^ 



^i+lj + ^i-l j ^ 



1 cos0 • (2d) 2 



-2(cos0j+ 2 

c °s0 j J cos0 j g 



.(£y ip+V\p- Vf) 



■] 



(23) 



with a relaxation tolerance of .1 meters. 
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One experiment, using the leapfrog scheme for time 
differencing, was performed with the "restorative-iterative" 
initialization method developed by V7inninghoff (1971) . 

This method involved using the Euler-backward time dif- 
ferencing scheme to alternately step forward and backward 
six times. After each iteration of equations (1), (2), 

and (3) the following restoration was added: 

(uh) i;j = (l-k u ) (uh) ±j + k u (uh) Q 

(vh) i j = (l-k v ) (vh) i j + k v (vh) Q (24) 

h ij = d'^h) Hi + HH 

where the k's are functions of latitude. k u and k v were 
.5 from latitude 20° S to 20° N, 0 from 40° N and S to the 
poles, and a linear variation between 0 and .5 between 
latitude 20° and 40°. k^ was .5 from 40° N and S to the 
poles, 0 between 20° S and 20° N, and a linear variation 
between 20° and 40°. 
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V. WAVE ANALYSIS METHOD 



To calculate the phase angles and amplitudes, a fourier 
analysis was performed at each five degrees of latitude 
around the latitude circle. The fourier series was ex- 
pressed as follows: 



F(x) = A Q + 53 (-km cos mx + s; *- n 

m 

= C o + 23 c m c °s (mx-<5 m ) 

m 



where 



B 



c m ~ 



m 



m 



sin(<$ m ) cos ( 6 m ) 



and 

_1 ®m 
6 m = tan 1 _ 

Since the input stream function involved only wave number 
six and a mean height, only C , Cg, and 6g values were ex- 
tracted from the fourier analysis. 
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VI . RESULTS 



All the experiments performed with the Arakawa averaging 
technique showed a considerable tilt backward at high lati- 
tudes in the phase propagation of the wave. This is to be 
expected since the smoothing of the gradients in the 
technique tends to slow down the rate of propagation. 

Gates (1959) has shown that, as the wavelength and Ax de- 
crease proportionally, the phase speed of the wave remains 
constant, and also that if the wavelength decreases and 
Ax remains constant, the phase speed will decrease relative 
to the exact value. In this model, the Arakawa averaging 
technique gives an effective Ax which is comparable to that 
at low latitudes, thus as the wavelength decreased toward 
the poles the phase speed also decreased. 

The result of this differential movement was to cause, 
eventually, the formation of closed highs and lows at the 
higher latitudes which propagated equatorward. It is 
believed that this instability is possible due to non- 
linear effects introduced after the field ceased to be 
harmonic in the latitudinal direction. 

The amplitudes in all the experiments showed a tendency 
to decrease at latitudes below 45° and increase above 45°. 
The mean height also tended to increase at the higher 
latitude (75° and above) . These amplitude variations are 
also believed to be caused by the nonlinear effects. All 
the methods, except the Euler-backward, had small amplitude 
gravity waves propagating with about a ten hour period. 
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Table I shows the comparison of the time required for 
a 120-hour forecast using each of the four methods. It 
also gives a comparison of the initial twenty-four phase 
speed for selected latitudes. 



TABLE I 



Time Differencing 
Method 


Time Required 
120hr Forecast 


Phase Speed for 


Initial 24hrs 


Equator 


30° 


60° | 


75° 


Leapforg 


32 


14.0 


10.7 


1.7 


-11.3 


Euler-Backward 


57 


13.3 


11.2 


1.5 


-9.7 


Leapfrog- 












Trapezoidal 


58 


13.7 


10.8 


1.5 


-11.5 


Adams-Bashford 


46 


14.3 


10.8 


1.8 


-11.2 



Note: Time in minutes 

Phase speed in degrees longitude per day 



The experiments performed without the averaging technique 
still showed a slight tendency to tilt backwards at the 
higher latitudes. This was not expected and was believed to 
have been caused by truncation errors due to special treat- 
ment near the poles using only a second-order difference 
approximation for the derivatives. This belief was based 
on the fact that the input stream function does not vary 
linearly near the poles, which caused problems earlier in 
the interpolation for (uh) and h . 

Table II gives a relative comparison of the time re- 
quired by each method using a 2.5-minute time step. 
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TABLE II 



Time Differencing 
Method 


Time Required in Minutes 


Leapfrog 


40 min. for a 32 hr. forecast 


Euler-Backward 


70 min. for a 30 hr. forecast 


Leapfrog-Trapezoidal 


70 min. for a 29 hr. forecast 


Adams-Bashford 


40 min. for a 32 hr. forecast 



A. RESULTS OF INDIVIDUAL EXPERIMENTS USING THE AVERAGING 

TECHNIQUE 

Experiment 1 . This experiment was performed using the 
leapfrog time differencing method. Figure 3 shows the phase 
angles as a function of latitude at twelve-hour intervals 
out to thirty-six hours and Fig. 8a shows phase angles at 
twenty-four hour intervals out to 120 hours. The amplitudes 
for wave number six and the mean heights are shown in Fig. 

13 for selected latitudes. 

Experiment 2 . The second experiment was performed using 
the leapfrog scheme and the "restorative- iterative" initial- 
ization method. This experiment was performed to see if the 
tilt of the phase lines could be reduced by letting the mass 
and wind fields "adjust" before performing the integrations. 
As can be seen in Fig. 4, the tilt was not reduced. 

Experiment 3 . The Euler-backward method was used for 
this experiment. The phase curves are shown in Figs. 5 and 
8b. The mean height and wave number six amplitudes are 
shown in Fig. 14. It should be noted that the gravity waves 
present in the other three methods are effectively damped 
out with this method. 



25 



Also the maximum variations in amplitudes, which are 
approximately equal for the other three methods , are 
slightly less since this scheme tends to damp all waves. 

' Experiment 4 . This experiment was performed using the 
leapfrog- trapezoidal method. Figures 6 and 8c show the 
phase angles vs latitude curves. The amplitudes are shown 
in Fig. 15. The largest gravity wave amplitudes were ob- 
served using this method. 

Experiment 5 . The last experiment using the Arakawa 
averaging scheme was performed using the Adams-Bashford 
method. Figures 7 and 8d show the phase relationships and 
Fig. 16 and amplitudes. There was very little difference 
in the results between this method and the leapfrog, except 
for the time required, see Table I, for the integrations. 

B. RESULTS OF INDIVIDUAL EXPERIMENTS WITHOUT THE AVERAGING 

TECHNIQUE 

Experiment 6 . This experiment, like the first experi- 
ment, was performed using the leapfrog scheme. The time 
step was reduced from 15 minutes to 2.5 minutes. The phase 
angle profiles are shown in Fig. 9 for zero, twelve, and 
twenty- four hours. 

Experiment 7 . This experiment was the same as experiment 
3, except At was 2.5 minutes. The phase angle results are 
shown in Fig. 10. 

Experiment 8 . This experiment was the same as experi- 
ment 4 with the exception of At, which was reduced to 2.5 
minutes. The phase profiles are shown in Fig. 11. 
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Experiment 9 . The same time differencing method was 
used as in experiment 5 . The time increment was reduced 
from 10 minutes to 2.5 minutes. Figure 12 shows the phase 
relationships for this experiment. 
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Latitude 



90 




FIGURE 3. PHASE ANGLE VS LATITUDE USING THE LEAPFROG 
SCHEME WITH THE ARAKAWA AVERAGING. 



Note: The input height field was constant at 85° lat. The 

phase angle at that latitude is the result of trunca- 
tion errors in the fourier analysis. 



28 




Latitude 




Longitude 



FIGURE 4. PHASE ANGLE VS LATITUDE USING THE LEAPFROG 
SCHEME WITH THE ARAKAWA AVERAGING AND 
WINNINGHOFF ' S "RESTORATIVE-ITERATIVE" 
INITIALIZATION. 
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.It 



Latitude 




FIGURE 5. PHASE ANGLE VS LATITUDE USING THE EULER- 

BACKWARD SCHEME WITH THE ARAKAWA AVERAGING. 

See note on Fig. 3. 
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Latitude 




FIGURE 6. PHASE ANGLE VS LATITUDE USING THE LEAPFROG- 

TRAPEZOIDAL SCHEME WITH THE ARAKAWA AVERAGING. 

See note on Fig. 3. 
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Latitude 




Longitude 

FIGURE 7. PHASE ANGLE VS LATITUDE USING THE ADAMS- 

BASHFORD SCHEME WITH THE ARAKAWA AVERAGING. 

See note on Fig. 3. 
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Latitude 



90 




SCHEME WITHOUT AVERAGING. 

See^note on Fig. 3. 

This 24 hour movement appears to be wrong compared to 
the averaging case but time did not permit a re-run to 
verify this movement. 



34 




FIGURE 10. PHASE ANGLE VS LATITUDE USING THE EULER- 
BACKWARD SCHEME WITHOUT AVERAGING. 



See note on Fig. 3. 
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FIGURE 11. PHASE ANGLE VS LATITUDE USING THE LEAPFROG- 
TRAPEZOIDAL SCHEME WITHOUT AVERAGING. 

See note on Fig. 3. 



36 



Latitude 




FIGURE 12. PHASE ANGLE VS LATITUDE USING THE ADAMS- 
BASHFORD SCHEME WITHOUT AVERAGING. 

See note on Fig. 3. 
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Time (Hours) 

1 Wave #6 Ampl + Mean Heigl 

Mean Height 

FIGURE 14. AMPLITUDE VS TIME FOR SELECTED LATITUDES USING THE EULER- BACKWARD 
SCHEME WITH AVERAGING. 
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FIGURE 15. AMPLITUDE VS TIME FOR SELECTED LATITUDES USING THE LEAPFROG- 
TRAPEZOIDAL SCHEME WITH AVERAGING. 
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VII. CONCLUSIONS 



The Arakawa averaging method caused some problems with 
the initial field which was represented by a spherical 
harmonic, but it is felt that with real data, where the 
longitudinal scale does not necessarily decrease with 
latitude, the method might not cause such severe problems. 
Considering the alternatives, such as a reduced time step, 
variable grid size, or variable time step, the Arawaka 
procedure is a simple and effective method for spherical 
prediction. The reduced time step is much too expensive 
in computer time to be practical. The abruptly changed 
grid size causes severe problems around the area of the 
change. A variable time step might prove to be acceptable 
but would involve some very complex programming. 

In the experiments performed, a second order one- 
dimensional interpolation was used since problems arose 
from using a two-dimensional linear interpolation and was 
the easiest to apply. In the real data cases, a two- 
dimensional second order interpolation scheme would prob- 
ably give better overall results . 

Overall the Euler-backward method gave the best results 
since it effectively reduced the amplitudes of the gravity 
waves, but was expensive in computer time. Considering 
time requirements and overall results, the leapfrog method 
is still the most desirable. Some further tests with com- 
binations of the methods might produce a method which gives 
good results and is acceptable as far as time required is 
concerned. 
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